Claisen Rearrangement with GFN1-xTB¶
Optimized transition state.¶
Requires: AMS2026 or later
Related documentation
This report summarizes a GFN1-xTB workflow for the Claisen rearrangement of C=CC(CCCC)C1C=CC(C)=CC1=O. The workflow consisted of reactant optimization with normal modes, a one-dimensional PES scan of the carbonyl oxygen to terminal vinyl carbon distance, a transition-state search started from the highest-energy scan point, and a product optimization from the shortest-distance scan structure.
Conclusions¶
Reaction free energy barrier from Gibbs free energies:
17.61 kcal/molReaction free energy from Gibbs free energies:
-16.51 kcal/molTransition-state verification:
1imaginary frequency/frequenciesImaginary frequency values (cm^-1):
-523.87
Scan Profile¶
The PES scan used the directly identified carbonyl oxygen and terminal vinyl carbon and sampled the O…C distance across 17 points, ending at 1.500 Angstrom.
Representative scan structures¶
Shortest-distance structure (1.500 Angstrom, -45.070213 Eh)
Longest-distance structure (4.736 Angstrom, -45.043347 Eh)
Highest-energy structure (1.904 Angstrom, -45.020033 Eh)
Optimized stationary points¶
Stationary Point |
Electronic Energy (Eh) |
Relative Electronic (kcal/mol) |
Gibbs Free Energy (Eh) |
Relative Gibbs (kcal/mol) |
|---|---|---|---|---|
Reactant |
-45.043346 |
0.000000 |
-44.789812 |
0.000000 |
TS |
-45.017765 |
16.052713 |
-44.761746 |
17.611167 |
Product |
-45.073091 |
-18.665161 |
-44.816126 |
-16.512827 |
Optimized reactant
Optimized transition state
Optimized product
Job provenance¶
reactant_opt: Optimized reactant with normal modes for the reference Gibbs free energy.pes_scan: Constrained PES scan used to generate the scan profile and TS starting geometry.ts_search: Transition-state search with initial Hessian and normal-mode verification.product_opt: Optimized product candidate taken from the shortest-distance scan structure.
Calculation inputs¶
reactant_opt¶
Properties
NormalModes yes
End
Task GeometryOptimization
System
Atoms
C 0.0352835768 -1.2634200451 2.0295182868
C -0.0415942574 -0.1145626957 1.3513671875
C 0.3768358030 -0.0107622019 -0.1023946895
C 1.3482601493 1.1778757204 -0.3037836994
C 2.7043903423 1.0005137599 0.4064036347
C 3.6037073244 -0.0421656093 -0.2748047891
C 4.9927169421 -0.0650315668 0.3569733474
C -0.8639092594 0.0835558950 -1.0353753505
C -1.6167867797 -1.2320759497 -1.0242134522
C -2.8481006207 -1.3596284649 -0.5174591373
C -3.5920200469 -0.2007888095 0.0203302547
C -4.9621855196 -0.4060955913 0.5962914008
C -3.0923854072 1.0353581165 -0.0982955957
C -1.7693748858 1.2549099789 -0.7024024375
O -1.3979590243 2.4318691017 -0.9545350183
H -0.2932245841 -1.3055311163 3.0627223262
H 0.4201319644 -2.1657625412 1.5646403918
H -0.4547906310 0.7594326185 1.8469683919
H 0.9219001025 -0.9336388195 -0.3976829884
H 0.8826924392 2.1060815765 0.0912568365
H 1.5283912179 1.3393815818 -1.3892591926
H 3.2253886735 1.9831444015 0.3897766333
H 2.5489641777 0.7270983739 1.4724197858
H 3.7064446671 0.1953577571 -1.3559394825
H 3.1616795189 -1.0552978497 -0.1752916690
H 4.9208180589 -0.3155000900 1.4369243466
H 5.6171362308 -0.8337967696 -0.1451860994
H 5.4865795359 0.9232818198 0.2413060231
H -0.5128916286 0.2263228659 -2.0809303623
H -1.1176209481 -2.1146362083 -1.4056763567
H -3.3119945930 -2.3382881694 -0.5095056928
H -5.0700115595 -1.4398065936 0.9877091215
H -5.7246495708 -0.2382537447 -0.1928669653
H -5.1447217319 0.2985644575 1.4359193845
H -3.6670996767 1.8922948115 0.2337513069
End
BondOrders
1 2 2.0
1 16 1.0
1 17 1.0
2 3 1.0
2 18 1.0
3 4 1.0
3 8 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.0
8 14 1.0
8 29 1.0
9 10 2.0
9 30 1.0
10 11 1.0
10 31 1.0
11 12 1.0
11 13 2.0
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.0
13 35 1.0
14 15 2.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
pes_scan¶
PESScan
ScanCoordinate
Distance 15 1 4.735989 1.500000
nPoints 17
End
End
Task PESScan
System
Atoms
C 0.0828356678 -1.2361609693 1.8600638635
C 0.0230169857 -0.1005602292 1.1840009238
C 0.3826447205 0.0205084136 -0.2724822300
C 1.3521611871 1.1876495589 -0.5166034737
C 2.6460120158 1.0536430092 0.2949181241
C 3.5226812171 -0.1192578462 -0.1586482087
C 4.8440952630 -0.1487866778 0.6087729338
C -0.9228364398 0.1942082855 -1.1223234263
C -1.6417212509 -1.1113136883 -1.1603288968
C -2.8348664870 -1.3002647840 -0.5845342490
C -3.5290877878 -0.2152724185 0.0761825981
C -4.8363082077 -0.5304791294 0.7410256425
C -3.0099669152 1.0319921635 0.0677756390
C -1.7445286569 1.3573294302 -0.5816422387
O -1.3421283198 2.4993510190 -0.6786401890
H -0.1867739182 -1.2949537692 2.9060301868
H 0.4054485891 -2.1640373293 1.4040032539
H -0.3132648021 0.8147181482 1.6634827407
H 0.8526673855 -0.9130807246 -0.6103291960
H 0.8652795062 2.1269807927 -0.2427291012
H 1.5917222736 1.2336264690 -1.5839435535
H 3.2135209300 1.9836685557 0.1837301315
H 2.4029556447 0.9387983613 1.3567719686
H 3.7270456881 -0.0291243509 -1.2313367049
H 2.9958873621 -1.0655744537 0.0032020358
H 4.6610192136 -0.2577301258 1.6789699085
H 5.4540419485 -0.9882669028 0.2733338621
H 5.4025750367 0.7741601198 0.4448678964
H -0.5977984126 0.4596345834 -2.1362258221
H -1.1214637929 -1.9293767174 -1.6486896873
H -3.3158183859 -2.2722781899 -0.6066649463
H -4.6890376443 -1.2952553745 1.5059415186
H -5.5423337805 -0.9250234583 0.0077423503
H -5.2667892044 0.3561142097 1.2052954354
H -3.5308866291 1.8644140190 0.5276865919
End
BondOrders
1 2 2.0
1 16 1.0
1 17 1.0
2 3 1.0
2 18 1.0
3 4 1.0
3 8 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.0
8 14 1.0
8 29 1.0
9 10 2.0
9 30 1.0
10 11 1.0
10 31 1.0
11 12 1.0
11 13 2.0
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.0
13 35 1.0
14 15 2.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
ts_search¶
GeometryOptimization
InitialHessian
Type Calculate
End
End
Properties
NormalModes yes
End
Task TransitionStateSearch
System
Atoms
C -0.1883855639 0.8849440844 2.8380677896
C 0.5627428136 1.3021798440 1.7761149239
C 0.2952333107 0.7218491339 0.4570370534
C 1.3447101401 1.0740177929 -0.5988077995
C 2.7398580096 0.5677441272 -0.2063178575
C 2.8468332670 -0.9605979608 -0.2429802578
C 4.2721905153 -1.4211851135 0.0591373277
C -1.1747522277 1.3440766657 -0.0038150205
C -1.5207378757 0.8630891391 -1.3508859600
C -2.3683829623 -0.1681694771 -1.5433814212
C -3.0696805006 -0.7696212742 -0.4470282606
C -3.9924790716 -1.9182697347 -0.7364163010
C -2.9158147501 -0.2700241146 0.8113791419
C -2.0845178261 0.8661991276 1.0809153734
O -1.9403895618 1.3267279730 2.2360511787
H -0.1493070037 1.4294235569 3.7730952063
H -0.4767901508 -0.1543179245 2.9358426632
H 0.9713072494 2.3092268737 1.7869621609
H 0.1389341018 -0.3610438632 0.4980558820
H 1.3877309740 2.1612653470 -0.7201423987
H 1.0669049571 0.6375254770 -1.5627741289
H 3.4679284037 0.9925994395 -0.9055578088
H 2.9913397275 0.9283248842 0.7965107426
H 2.5493453999 -1.3238854757 -1.2328140760
H 2.1656346027 -1.3999529135 0.4928224177
H 4.5804954685 -1.0829907202 1.0494557149
H 4.3281249818 -2.5097136723 0.0308707210
H 4.9691052791 -1.0181568941 -0.6771753844
H -1.0207148970 2.4248273894 0.0189637741
H -1.0450925835 1.3493498217 -2.1957393948
H -2.5748585113 -0.5292833208 -2.5455309316
H -3.4370217067 -2.7390101685 -1.1940103991
H -4.7684533300 -1.6087508459 -1.4392993575
H -4.4676411820 -2.2800376390 0.1748340102
H -3.4833994970 -0.6683595645 1.6452363582
End
BondOrders
1 2 2.0
1 16 1.0
1 17 1.0
2 3 1.0
2 18 1.0
3 4 1.0
3 8 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.0
8 14 1.0
8 29 1.0
9 10 2.0
9 30 1.0
10 11 1.0
10 31 1.0
11 12 1.0
11 13 2.0
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.0
13 35 1.0
14 15 2.0
End
End
System final
Atoms
C -1.0451269478 2.0735889566 2.2034291742
C 0.3809633753 1.9201215848 1.7721373932
C 1.1744665724 0.9376768117 2.1784859003
C 2.5918951088 0.7657262339 1.7151237820
C 2.7112172482 -0.4246946392 0.7469721264
C 1.9541349796 -0.1700007854 -0.5601651801
C 2.0736467518 -1.3557249519 -1.5155867636
C -2.6204936975 1.4665200110 -1.1058903775
C -2.7938208214 0.4675898504 -2.0476039400
C -2.3983555033 -0.8358684068 -1.7758915240
C -1.8296740847 -1.1465175506 -0.5444949898
C -1.4268441422 -2.5633638543 -0.2308997241
C -1.6530482167 -0.1460785668 0.4082593299
C -2.0439493279 1.1630869836 0.1305154182
O -1.9305596711 2.1983235634 0.9990828078
H -1.2003629359 3.0217346306 2.7262971868
H -1.3943747537 1.2462085997 2.8296833104
H 0.7319447437 2.6783949005 1.0785449878
H 0.8086725407 0.1766513728 2.8640131241
H 3.2361043614 0.5779133975 2.5800608176
H 2.9429714280 1.6751678780 1.2168056076
H 3.7705997359 -0.5934415113 0.5290325792
H 2.3197285611 -1.3262388519 1.2303528490
H 2.3530861556 0.7287477483 -1.0429599793
H 0.8970543304 0.0117095184 -0.3400784601
H 1.6755236766 -2.2602423039 -1.0532884067
H 1.5119163852 -1.1570786327 -2.4287278953
H 3.1172183154 -1.5327563377 -1.7801746401
H -2.9237811573 2.4862797052 -1.3082386669
H -3.2380402552 0.7066147322 -3.0077839698
H -2.5339668703 -1.6104563479 -2.5222966220
H -1.0206620652 -3.0537871017 -1.1160991236
H -2.2983720908 -3.1326597682 0.1005183787
H -0.6800413687 -2.5911646697 0.5628480847
H -1.2196703602 -0.4019821983 1.3666930871
End
BondOrders
1 2 1.0
1 15 1.0
1 16 1.0
1 17 1.0
2 3 2.0
2 18 1.0
3 4 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.5
8 14 1.5
8 29 1.0
9 10 1.5
9 30 1.0
10 11 1.5
10 31 1.0
11 12 1.0
11 13 1.5
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.5
13 35 1.0
14 15 1.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
product_opt¶
Properties
NormalModes yes
End
Task GeometryOptimization
System
Atoms
C -1.0451269478 2.0735889566 2.2034291742
C 0.3809633753 1.9201215848 1.7721373932
C 1.1744665724 0.9376768117 2.1784859003
C 2.5918951088 0.7657262339 1.7151237820
C 2.7112172482 -0.4246946392 0.7469721264
C 1.9541349796 -0.1700007854 -0.5601651801
C 2.0736467518 -1.3557249519 -1.5155867636
C -2.6204936975 1.4665200110 -1.1058903775
C -2.7938208214 0.4675898504 -2.0476039400
C -2.3983555033 -0.8358684068 -1.7758915240
C -1.8296740847 -1.1465175506 -0.5444949898
C -1.4268441422 -2.5633638543 -0.2308997241
C -1.6530482167 -0.1460785668 0.4082593299
C -2.0439493279 1.1630869836 0.1305154182
O -1.9305596711 2.1983235634 0.9990828078
H -1.2003629359 3.0217346306 2.7262971868
H -1.3943747537 1.2462085997 2.8296833104
H 0.7319447437 2.6783949005 1.0785449878
H 0.8086725407 0.1766513728 2.8640131241
H 3.2361043614 0.5779133975 2.5800608176
H 2.9429714280 1.6751678780 1.2168056076
H 3.7705997359 -0.5934415113 0.5290325792
H 2.3197285611 -1.3262388519 1.2303528490
H 2.3530861556 0.7287477483 -1.0429599793
H 0.8970543304 0.0117095184 -0.3400784601
H 1.6755236766 -2.2602423039 -1.0532884067
H 1.5119163852 -1.1570786327 -2.4287278953
H 3.1172183154 -1.5327563377 -1.7801746401
H -2.9237811573 2.4862797052 -1.3082386669
H -3.2380402552 0.7066147322 -3.0077839698
H -2.5339668703 -1.6104563479 -2.5222966220
H -1.0206620652 -3.0537871017 -1.1160991236
H -2.2983720908 -3.1326597682 0.1005183787
H -0.6800413687 -2.5911646697 0.5628480847
H -1.2196703602 -0.4019821983 1.3666930871
End
BondOrders
1 2 1.0
1 15 1.0
1 16 1.0
1 17 1.0
2 3 2.0
2 18 1.0
3 4 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.5
8 14 1.5
8 29 1.0
9 10 1.5
9 30 1.0
10 11 1.5
10 31 1.0
11 12 1.0
11 13 1.5
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.5
13 35 1.0
14 15 1.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
Prompts and Python scripts¶
Prompt (instruction for AI agent)
Use $ams2026
I want the reaction free energy barrier and reaction free energies for the
Claisen rearrangement of the molecule with SMILES C=CC(CCCC)C1C=CC(C)=CC1=O
Find the carbonyl oxygen and terminal vinyl carbon ChemicalSystem bonds methods
Method: GFN1-xTB
Do a bond scan where that distance is scanned to 1.5 angstrom in intervals of
0.2 angstrom
Plot the energy vs distance (include in report), including an images of the
structures with shortest distance, longest distance, and the highest energy.
Use the highest energy point as initial structure for transition state search.
Make sure to calculate the initial hessian for the TS search and the normal
modes for free energy. Verify that TS has exactly one imaginary (negative)
frequency. Also optimize the reactant and product with normal modes.
Include in the report pictures and energies and relative energies (in kcal/mol)
of optimized reactant, TS, and product.
01-run.py
#!/usr/bin/env amspython
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
import scm.plams as plams
from scm.base import ChemicalSystem
from scm.input_classes import AMS
from scm.plams import AMSJob, Settings
SMILES = "C=CC(CCCC)C1C=CC(C)=CC1=O"
TARGET_DISTANCE_ANGSTROM = 1.5
SCAN_STEP_ANGSTROM = 0.2
@dataclass(frozen=True)
class AtomPair:
carbonyl_oxygen: int
terminal_vinyl_carbon: int
def ensure_bonds(system: object) -> None:
has_bonds = getattr(system, "has_bonds", None)
if callable(has_bonds):
if not has_bonds():
system.guess_bonds()
return
bonds = getattr(system, "bonds", None)
if bonds is not None and len(bonds) == 0:
system.guess_bonds()
def atom_symbol(system: ChemicalSystem, index: int) -> str:
return system.atoms[index].element.symbol
def adjacency(system: ChemicalSystem) -> dict[int, list[tuple[int, float]]]:
graph: dict[int, list[tuple[int, float]]] = {i: [] for i in range(len(system.atoms))}
for i, j, bond in system.bonds:
graph[i].append((j, bond.order))
graph[j].append((i, bond.order))
return graph
def heavy_neighbors(system: ChemicalSystem, graph: dict[int, list[tuple[int, float]]], index: int) -> list[tuple[int, float]]:
return [(nbr, order) for nbr, order in graph[index] if atom_symbol(system, nbr) != "H"]
def distance_angstrom(system: ChemicalSystem, i: int, j: int) -> float:
coords_i = np.array(system.atoms[i].coords)
coords_j = np.array(system.atoms[j].coords)
return float(np.linalg.norm(coords_i - coords_j))
def find_target_atoms(system: ChemicalSystem) -> AtomPair:
ensure_bonds(system)
graph = adjacency(system)
carbonyl_oxygen_candidates: list[int] = []
for index in range(len(system.atoms)):
if atom_symbol(system, index) != "O":
continue
nbrs = heavy_neighbors(system, graph, index)
if len(nbrs) != 1:
continue
nbr, order = nbrs[0]
if atom_symbol(system, nbr) == "C" and order >= 1.5:
carbonyl_oxygen_candidates.append(index)
if len(carbonyl_oxygen_candidates) != 1:
raise RuntimeError(f"Expected exactly one carbonyl oxygen, found {carbonyl_oxygen_candidates}")
carbonyl_oxygen = carbonyl_oxygen_candidates[0]
terminal_vinyl_candidates: list[int] = []
for index in range(len(system.atoms)):
if atom_symbol(system, index) != "C":
continue
nbrs = heavy_neighbors(system, graph, index)
if len(nbrs) != 1:
continue
nbr, order = nbrs[0]
if atom_symbol(system, nbr) == "C" and order >= 1.5:
terminal_vinyl_candidates.append(index)
if not terminal_vinyl_candidates:
raise RuntimeError("Could not find a terminal vinyl carbon from the bond graph")
terminal_vinyl_carbon = max(
terminal_vinyl_candidates,
key=lambda idx: distance_angstrom(system, carbonyl_oxygen, idx),
)
return AtomPair(carbonyl_oxygen=carbonyl_oxygen, terminal_vinyl_carbon=terminal_vinyl_carbon)
def validate(settings: Settings) -> None:
AMS.from_settings(settings)
def dftb_settings() -> Settings:
settings = Settings()
settings.input.DFTB.Model = "GFN1-xTB"
return settings
def reactant_settings() -> Settings:
settings = dftb_settings()
settings.input.ams.Task = "GeometryOptimization"
settings.input.ams.Properties.NormalModes = "Yes"
validate(settings)
return settings
def pesscan_settings(atom_pair: AtomPair, start_distance: float, end_distance: float, npoints: int) -> Settings:
settings = dftb_settings()
settings.input.ams.Task = "PESScan"
settings.input.ams.PESScan.ScanCoordinate.nPoints = npoints
settings.input.ams.PESScan.ScanCoordinate.Distance = (
f"{atom_pair.carbonyl_oxygen + 1:d} {atom_pair.terminal_vinyl_carbon + 1:d} "
f"{start_distance:.6f} {end_distance:.6f}"
)
validate(settings)
return settings
def ts_settings() -> Settings:
settings = dftb_settings()
settings.input.ams.Task = "TransitionStateSearch"
settings.input.ams.Properties.NormalModes = "Yes"
settings.input.ams.GeometryOptimization.InitialHessian.Type = "Calculate"
validate(settings)
return settings
def product_settings() -> Settings:
settings = dftb_settings()
settings.input.ams.Task = "GeometryOptimization"
settings.input.ams.Properties.NormalModes = "Yes"
validate(settings)
return settings
def run_job(name: str, molecule: ChemicalSystem | dict[str, ChemicalSystem], settings: Settings) -> AMSJob:
job = AMSJob(molecule=molecule, settings=settings, name=name)
result = job.run()
if not result.ok():
raise RuntimeError(f"{name} failed")
return job
def main() -> None:
plams.init(folder="01-run_workdir")
system = ChemicalSystem.from_smiles(SMILES)
ensure_bonds(system)
initial_atom_pair = find_target_atoms(system)
print(
"Initial atom indices (1-based):",
f"carbonyl O={initial_atom_pair.carbonyl_oxygen + 1}",
f"terminal vinyl C={initial_atom_pair.terminal_vinyl_carbon + 1}",
)
reactant_job = run_job("reactant_opt", system, reactant_settings())
reactant_system = reactant_job.results.get_main_system()
ensure_bonds(reactant_system)
atom_pair = find_target_atoms(reactant_system)
start_distance = distance_angstrom(reactant_system, atom_pair.carbonyl_oxygen, atom_pair.terminal_vinyl_carbon)
if start_distance <= TARGET_DISTANCE_ANGSTROM:
raise RuntimeError(
f"Reactant O...C distance is already {start_distance:.3f} A, cannot scan down to {TARGET_DISTANCE_ANGSTROM:.3f} A"
)
npoints = int(round((start_distance - TARGET_DISTANCE_ANGSTROM) / SCAN_STEP_ANGSTROM)) + 1
if npoints < 2:
raise RuntimeError(f"Invalid PES scan point count from start distance {start_distance:.3f} A")
pesscan_job = run_job(
"pes_scan",
reactant_system,
pesscan_settings(atom_pair, start_distance, TARGET_DISTANCE_ANGSTROM, npoints),
)
pesscan_results = pesscan_job.results.get_pesscan_results()
energies = np.array(pesscan_results["PES"], dtype=float)
distances = np.array(pesscan_results["RaveledPESCoords"][0], dtype=float)
molecules = pesscan_results["Molecules"]
highest_index = int(np.argmax(energies))
shortest_index = int(np.argmin(distances))
highest_system = molecules[highest_index]
product_guess = molecules[shortest_index]
ensure_bonds(highest_system)
ensure_bonds(product_guess)
ts_job = run_job(
"ts_search",
{"": highest_system, "final": product_guess},
ts_settings(),
)
ts_frequencies = np.array(ts_job.results.get_frequencies(unit="cm^-1"), dtype=float)
n_imaginary = int(np.count_nonzero(ts_frequencies < 0.0))
print(f"Transition state imaginary frequencies: {n_imaginary}")
if n_imaginary != 1:
raise RuntimeError(f"Transition state verification failed: expected 1 imaginary frequency, found {n_imaginary}")
product_job = run_job("product_opt", product_guess, product_settings())
product_system = product_job.results.get_main_system()
ensure_bonds(product_system)
print("Finished workflow")
print(
f"Final atom indices (1-based) used for scan: O={atom_pair.carbonyl_oxygen + 1}, "
f"C={atom_pair.terminal_vinyl_carbon + 1}"
)
print(f"Scan start distance: {start_distance:.3f} A")
print(f"Requested scan end distance: {TARGET_DISTANCE_ANGSTROM:.3f} A")
print(f"Number of PES points: {npoints}")
if __name__ == "__main__":
main()
report.py
#!/usr/bin/env amspython
from __future__ import annotations
from dataclasses import dataclass
from pathlib import Path
import shutil
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scm.plams as plams
from scm.base import ChemicalSystem, Units
from scm.plams import AMSJob
SMILES = "C=CC(CCCC)C1C=CC(C)=CC1=O"
WORKDIR = Path("01-run_workdir")
REPORT_PATH = Path("report.md")
BACKUP_PATH = Path("report.md.bk")
@dataclass(frozen=True)
class LoadedJob:
job: AMSJob
why: str
def load_jobs() -> dict[str, LoadedJob]:
specs = {
"reactant": LoadedJob(
job=AMSJob.load_external(str(WORKDIR / "reactant_opt")),
why="Optimized reactant with normal modes for the reference Gibbs free energy.",
),
"scan": LoadedJob(
job=AMSJob.load_external(str(WORKDIR / "pes_scan")),
why="Constrained PES scan used to generate the scan profile and TS starting geometry.",
),
"ts": LoadedJob(
job=AMSJob.load_external(str(WORKDIR / "ts_search")),
why="Transition-state search with initial Hessian and normal-mode verification.",
),
"product": LoadedJob(
job=AMSJob.load_external(str(WORKDIR / "product_opt")),
why="Optimized product candidate taken from the shortest-distance scan structure.",
),
}
return specs
def energy_hartree(job: AMSJob) -> float:
return float(job.results.get_energy(unit="hartree"))
def gibbs_hartree(job: AMSJob) -> float:
return float(job.results.readrkf("Thermodynamics", "Gibbs free Energy", file="engine"))
def relative_kcal(values_hartree: dict[str, float], reference_key: str) -> dict[str, float]:
factor = Units.conversion_factor("hartree", "kcal/mol")
reference = values_hartree[reference_key]
return {key: (value - reference) * factor for key, value in values_hartree.items()}
def render_system_image(system: ChemicalSystem, picture_path: Path) -> None:
plams.view(
system,
guess_bonds=True,
direction="tilt_pca3",
width=360,
height=280,
picture_path=str(picture_path),
)
def write_report(markdown: str) -> None:
if REPORT_PATH.exists():
shutil.copyfile(REPORT_PATH, BACKUP_PATH)
REPORT_PATH.write_text(markdown, encoding="utf-8")
def main() -> None:
jobs = load_jobs()
reactant_job = jobs["reactant"].job
scan_job = jobs["scan"].job
ts_job = jobs["ts"].job
product_job = jobs["product"].job
scan_results = scan_job.results.get_pesscan_results()
scan_energies = np.array(scan_results["PES"], dtype=float)
scan_distances_bohr = np.array(scan_results["RaveledPESCoords"][0], dtype=float)
scan_distances = scan_distances_bohr * Units.conversion_factor("bohr", "angstrom")
scan_molecules = scan_results["Molecules"]
highest_index = int(np.argmax(scan_energies))
shortest_index = int(np.argmin(scan_distances))
longest_index = int(np.argmax(scan_distances))
energy_plot_path = Path("scan_energy_vs_distance.png")
shortest_image_path = Path("scan_shortest_distance.png")
longest_image_path = Path("scan_longest_distance.png")
highest_image_path = Path("scan_highest_energy.png")
reactant_image_path = Path("reactant_optimized.png")
ts_image_path = Path("ts_optimized.png")
product_image_path = Path("product_optimized.png")
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(scan_distances, scan_energies, marker="o")
ax.set_xlabel("O...C distance (Angstrom)")
ax.set_ylabel("Energy (hartree)")
ax.set_title("GFN1-xTB PES scan")
fig.tight_layout()
fig.savefig(energy_plot_path, dpi=200)
plt.close(fig)
render_system_image(scan_molecules[shortest_index], shortest_image_path)
render_system_image(scan_molecules[longest_index], longest_image_path)
render_system_image(scan_molecules[highest_index], highest_image_path)
render_system_image(reactant_job.results.get_main_system(), reactant_image_path)
render_system_image(ts_job.results.get_main_system(), ts_image_path)
render_system_image(product_job.results.get_main_system(), product_image_path)
electronic = {
"Reactant": energy_hartree(reactant_job),
"TS": energy_hartree(ts_job),
"Product": energy_hartree(product_job),
}
gibbs = {
"Reactant": gibbs_hartree(reactant_job),
"TS": gibbs_hartree(ts_job),
"Product": gibbs_hartree(product_job),
}
rel_electronic = relative_kcal(electronic, "Reactant")
rel_gibbs = relative_kcal(gibbs, "Reactant")
summary = pd.DataFrame(
{
"Electronic Energy (Eh)": electronic,
"Relative Electronic (kcal/mol)": rel_electronic,
"Gibbs Free Energy (Eh)": gibbs,
"Relative Gibbs (kcal/mol)": rel_gibbs,
}
)
summary.index.name = "Stationary Point"
ts_frequencies = np.array(ts_job.results.get_frequencies(unit="cm^-1"), dtype=float)
imaginary_frequencies = ts_frequencies[ts_frequencies < 0.0]
barrier_g_kcal = rel_gibbs["TS"]
reaction_free_energy_kcal = rel_gibbs["Product"]
markdown = f"""# Claisen Rearrangement with GFN1-xTB
This report summarizes a GFN1-xTB workflow for the Claisen rearrangement of `{SMILES}`. The workflow consisted of reactant optimization with normal modes, a one-dimensional PES scan of the carbonyl oxygen to terminal vinyl carbon distance, a transition-state search started from the highest-energy scan point, and a product optimization from the shortest-distance scan structure.
## Conclusions
- Reaction free energy barrier from Gibbs free energies: `{barrier_g_kcal:.2f} kcal/mol`
- Reaction free energy from Gibbs free energies: `{reaction_free_energy_kcal:.2f} kcal/mol`
- Transition-state verification: `{len(imaginary_frequencies)}` imaginary frequency/frequencies
- Imaginary frequency values (cm^-1): `{", ".join(f"{value:.2f}" for value in imaginary_frequencies)}`
## Scan Profile
The PES scan used the directly identified carbonyl oxygen and terminal vinyl carbon and sampled the O...C distance across `{len(scan_distances)}` points, ending at `{scan_distances[shortest_index]:.3f} Angstrom`.
### Representative scan structures
Shortest-distance structure (`{scan_distances[shortest_index]:.3f} Angstrom`, `{scan_energies[shortest_index]:.6f} Eh`)
Longest-distance structure (`{scan_distances[longest_index]:.3f} Angstrom`, `{scan_energies[longest_index]:.6f} Eh`)
Highest-energy structure (`{scan_distances[highest_index]:.3f} Angstrom`, `{scan_energies[highest_index]:.6f} Eh`)
## Optimized stationary points
{summary.to_markdown(floatfmt=".6f")}
Optimized reactant
Optimized transition state
Optimized product
## Job provenance
- `reactant_opt`: {jobs["reactant"].why}
- `pes_scan`: {jobs["scan"].why}
- `ts_search`: {jobs["ts"].why}
- `product_opt`: {jobs["product"].why}
## Calculation inputs
### `reactant_opt`
```text
{reactant_job.get_input()}
```
### `pes_scan`
```text
{scan_job.get_input()}
```
### `ts_search`
```text
{ts_job.get_input()}
```
### `product_opt`
```text
{product_job.get_input()}
```
"""
write_report(markdown)
if __name__ == "__main__":
main()
Original Markdown report
# Claisen Rearrangement with GFN1-xTB
This report summarizes a GFN1-xTB workflow for the Claisen rearrangement of `C=CC(CCCC)C1C=CC(C)=CC1=O`. The workflow consisted of reactant optimization with normal modes, a one-dimensional PES scan of the carbonyl oxygen to terminal vinyl carbon distance, a transition-state search started from the highest-energy scan point, and a product optimization from the shortest-distance scan structure.
## Conclusions
- Reaction free energy barrier from Gibbs free energies: `17.61 kcal/mol`
- Reaction free energy from Gibbs free energies: `-16.51 kcal/mol`
- Transition-state verification: `1` imaginary frequency/frequencies
- Imaginary frequency values (cm^-1): `-523.87`
## Scan Profile
The PES scan used the directly identified carbonyl oxygen and terminal vinyl carbon and sampled the O...C distance across `17` points, ending at `1.500 Angstrom`.
### Representative scan structures
Shortest-distance structure (`1.500 Angstrom`, `-45.070213 Eh`)
Longest-distance structure (`4.736 Angstrom`, `-45.043347 Eh`)
Highest-energy structure (`1.904 Angstrom`, `-45.020033 Eh`)
## Optimized stationary points
| Stationary Point | Electronic Energy (Eh) | Relative Electronic (kcal/mol) | Gibbs Free Energy (Eh) | Relative Gibbs (kcal/mol) |
|:-------------------|-------------------------:|---------------------------------:|-------------------------:|----------------------------:|
| Reactant | -45.043346 | 0.000000 | -44.789812 | 0.000000 |
| TS | -45.017765 | 16.052713 | -44.761746 | 17.611167 |
| Product | -45.073091 | -18.665161 | -44.816126 | -16.512827 |
Optimized reactant
Optimized transition state
Optimized product
## Job provenance
- `reactant_opt`: Optimized reactant with normal modes for the reference Gibbs free energy.
- `pes_scan`: Constrained PES scan used to generate the scan profile and TS starting geometry.
- `ts_search`: Transition-state search with initial Hessian and normal-mode verification.
- `product_opt`: Optimized product candidate taken from the shortest-distance scan structure.
## Calculation inputs
### `reactant_opt`
```ams
Properties
NormalModes yes
End
Task GeometryOptimization
System
Atoms
C 0.0352835768 -1.2634200451 2.0295182868
C -0.0415942574 -0.1145626957 1.3513671875
C 0.3768358030 -0.0107622019 -0.1023946895
C 1.3482601493 1.1778757204 -0.3037836994
C 2.7043903423 1.0005137599 0.4064036347
C 3.6037073244 -0.0421656093 -0.2748047891
C 4.9927169421 -0.0650315668 0.3569733474
C -0.8639092594 0.0835558950 -1.0353753505
C -1.6167867797 -1.2320759497 -1.0242134522
C -2.8481006207 -1.3596284649 -0.5174591373
C -3.5920200469 -0.2007888095 0.0203302547
C -4.9621855196 -0.4060955913 0.5962914008
C -3.0923854072 1.0353581165 -0.0982955957
C -1.7693748858 1.2549099789 -0.7024024375
O -1.3979590243 2.4318691017 -0.9545350183
H -0.2932245841 -1.3055311163 3.0627223262
H 0.4201319644 -2.1657625412 1.5646403918
H -0.4547906310 0.7594326185 1.8469683919
H 0.9219001025 -0.9336388195 -0.3976829884
H 0.8826924392 2.1060815765 0.0912568365
H 1.5283912179 1.3393815818 -1.3892591926
H 3.2253886735 1.9831444015 0.3897766333
H 2.5489641777 0.7270983739 1.4724197858
H 3.7064446671 0.1953577571 -1.3559394825
H 3.1616795189 -1.0552978497 -0.1752916690
H 4.9208180589 -0.3155000900 1.4369243466
H 5.6171362308 -0.8337967696 -0.1451860994
H 5.4865795359 0.9232818198 0.2413060231
H -0.5128916286 0.2263228659 -2.0809303623
H -1.1176209481 -2.1146362083 -1.4056763567
H -3.3119945930 -2.3382881694 -0.5095056928
H -5.0700115595 -1.4398065936 0.9877091215
H -5.7246495708 -0.2382537447 -0.1928669653
H -5.1447217319 0.2985644575 1.4359193845
H -3.6670996767 1.8922948115 0.2337513069
End
BondOrders
1 2 2.0
1 16 1.0
1 17 1.0
2 3 1.0
2 18 1.0
3 4 1.0
3 8 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.0
8 14 1.0
8 29 1.0
9 10 2.0
9 30 1.0
10 11 1.0
10 31 1.0
11 12 1.0
11 13 2.0
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.0
13 35 1.0
14 15 2.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
```
### `pes_scan`
```ams
PESScan
ScanCoordinate
Distance 15 1 4.735989 1.500000
nPoints 17
End
End
Task PESScan
System
Atoms
C 0.0828356678 -1.2361609693 1.8600638635
C 0.0230169857 -0.1005602292 1.1840009238
C 0.3826447205 0.0205084136 -0.2724822300
C 1.3521611871 1.1876495589 -0.5166034737
C 2.6460120158 1.0536430092 0.2949181241
C 3.5226812171 -0.1192578462 -0.1586482087
C 4.8440952630 -0.1487866778 0.6087729338
C -0.9228364398 0.1942082855 -1.1223234263
C -1.6417212509 -1.1113136883 -1.1603288968
C -2.8348664870 -1.3002647840 -0.5845342490
C -3.5290877878 -0.2152724185 0.0761825981
C -4.8363082077 -0.5304791294 0.7410256425
C -3.0099669152 1.0319921635 0.0677756390
C -1.7445286569 1.3573294302 -0.5816422387
O -1.3421283198 2.4993510190 -0.6786401890
H -0.1867739182 -1.2949537692 2.9060301868
H 0.4054485891 -2.1640373293 1.4040032539
H -0.3132648021 0.8147181482 1.6634827407
H 0.8526673855 -0.9130807246 -0.6103291960
H 0.8652795062 2.1269807927 -0.2427291012
H 1.5917222736 1.2336264690 -1.5839435535
H 3.2135209300 1.9836685557 0.1837301315
H 2.4029556447 0.9387983613 1.3567719686
H 3.7270456881 -0.0291243509 -1.2313367049
H 2.9958873621 -1.0655744537 0.0032020358
H 4.6610192136 -0.2577301258 1.6789699085
H 5.4540419485 -0.9882669028 0.2733338621
H 5.4025750367 0.7741601198 0.4448678964
H -0.5977984126 0.4596345834 -2.1362258221
H -1.1214637929 -1.9293767174 -1.6486896873
H -3.3158183859 -2.2722781899 -0.6066649463
H -4.6890376443 -1.2952553745 1.5059415186
H -5.5423337805 -0.9250234583 0.0077423503
H -5.2667892044 0.3561142097 1.2052954354
H -3.5308866291 1.8644140190 0.5276865919
End
BondOrders
1 2 2.0
1 16 1.0
1 17 1.0
2 3 1.0
2 18 1.0
3 4 1.0
3 8 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.0
8 14 1.0
8 29 1.0
9 10 2.0
9 30 1.0
10 11 1.0
10 31 1.0
11 12 1.0
11 13 2.0
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.0
13 35 1.0
14 15 2.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
```
### `ts_search`
```ams
GeometryOptimization
InitialHessian
Type Calculate
End
End
Properties
NormalModes yes
End
Task TransitionStateSearch
System
Atoms
C -0.1883855639 0.8849440844 2.8380677896
C 0.5627428136 1.3021798440 1.7761149239
C 0.2952333107 0.7218491339 0.4570370534
C 1.3447101401 1.0740177929 -0.5988077995
C 2.7398580096 0.5677441272 -0.2063178575
C 2.8468332670 -0.9605979608 -0.2429802578
C 4.2721905153 -1.4211851135 0.0591373277
C -1.1747522277 1.3440766657 -0.0038150205
C -1.5207378757 0.8630891391 -1.3508859600
C -2.3683829623 -0.1681694771 -1.5433814212
C -3.0696805006 -0.7696212742 -0.4470282606
C -3.9924790716 -1.9182697347 -0.7364163010
C -2.9158147501 -0.2700241146 0.8113791419
C -2.0845178261 0.8661991276 1.0809153734
O -1.9403895618 1.3267279730 2.2360511787
H -0.1493070037 1.4294235569 3.7730952063
H -0.4767901508 -0.1543179245 2.9358426632
H 0.9713072494 2.3092268737 1.7869621609
H 0.1389341018 -0.3610438632 0.4980558820
H 1.3877309740 2.1612653470 -0.7201423987
H 1.0669049571 0.6375254770 -1.5627741289
H 3.4679284037 0.9925994395 -0.9055578088
H 2.9913397275 0.9283248842 0.7965107426
H 2.5493453999 -1.3238854757 -1.2328140760
H 2.1656346027 -1.3999529135 0.4928224177
H 4.5804954685 -1.0829907202 1.0494557149
H 4.3281249818 -2.5097136723 0.0308707210
H 4.9691052791 -1.0181568941 -0.6771753844
H -1.0207148970 2.4248273894 0.0189637741
H -1.0450925835 1.3493498217 -2.1957393948
H -2.5748585113 -0.5292833208 -2.5455309316
H -3.4370217067 -2.7390101685 -1.1940103991
H -4.7684533300 -1.6087508459 -1.4392993575
H -4.4676411820 -2.2800376390 0.1748340102
H -3.4833994970 -0.6683595645 1.6452363582
End
BondOrders
1 2 2.0
1 16 1.0
1 17 1.0
2 3 1.0
2 18 1.0
3 4 1.0
3 8 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.0
8 14 1.0
8 29 1.0
9 10 2.0
9 30 1.0
10 11 1.0
10 31 1.0
11 12 1.0
11 13 2.0
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.0
13 35 1.0
14 15 2.0
End
End
System final
Atoms
C -1.0451269478 2.0735889566 2.2034291742
C 0.3809633753 1.9201215848 1.7721373932
C 1.1744665724 0.9376768117 2.1784859003
C 2.5918951088 0.7657262339 1.7151237820
C 2.7112172482 -0.4246946392 0.7469721264
C 1.9541349796 -0.1700007854 -0.5601651801
C 2.0736467518 -1.3557249519 -1.5155867636
C -2.6204936975 1.4665200110 -1.1058903775
C -2.7938208214 0.4675898504 -2.0476039400
C -2.3983555033 -0.8358684068 -1.7758915240
C -1.8296740847 -1.1465175506 -0.5444949898
C -1.4268441422 -2.5633638543 -0.2308997241
C -1.6530482167 -0.1460785668 0.4082593299
C -2.0439493279 1.1630869836 0.1305154182
O -1.9305596711 2.1983235634 0.9990828078
H -1.2003629359 3.0217346306 2.7262971868
H -1.3943747537 1.2462085997 2.8296833104
H 0.7319447437 2.6783949005 1.0785449878
H 0.8086725407 0.1766513728 2.8640131241
H 3.2361043614 0.5779133975 2.5800608176
H 2.9429714280 1.6751678780 1.2168056076
H 3.7705997359 -0.5934415113 0.5290325792
H 2.3197285611 -1.3262388519 1.2303528490
H 2.3530861556 0.7287477483 -1.0429599793
H 0.8970543304 0.0117095184 -0.3400784601
H 1.6755236766 -2.2602423039 -1.0532884067
H 1.5119163852 -1.1570786327 -2.4287278953
H 3.1172183154 -1.5327563377 -1.7801746401
H -2.9237811573 2.4862797052 -1.3082386669
H -3.2380402552 0.7066147322 -3.0077839698
H -2.5339668703 -1.6104563479 -2.5222966220
H -1.0206620652 -3.0537871017 -1.1160991236
H -2.2983720908 -3.1326597682 0.1005183787
H -0.6800413687 -2.5911646697 0.5628480847
H -1.2196703602 -0.4019821983 1.3666930871
End
BondOrders
1 2 1.0
1 15 1.0
1 16 1.0
1 17 1.0
2 3 2.0
2 18 1.0
3 4 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.5
8 14 1.5
8 29 1.0
9 10 1.5
9 30 1.0
10 11 1.5
10 31 1.0
11 12 1.0
11 13 1.5
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.5
13 35 1.0
14 15 1.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
```
### `product_opt`
```ams
Properties
NormalModes yes
End
Task GeometryOptimization
System
Atoms
C -1.0451269478 2.0735889566 2.2034291742
C 0.3809633753 1.9201215848 1.7721373932
C 1.1744665724 0.9376768117 2.1784859003
C 2.5918951088 0.7657262339 1.7151237820
C 2.7112172482 -0.4246946392 0.7469721264
C 1.9541349796 -0.1700007854 -0.5601651801
C 2.0736467518 -1.3557249519 -1.5155867636
C -2.6204936975 1.4665200110 -1.1058903775
C -2.7938208214 0.4675898504 -2.0476039400
C -2.3983555033 -0.8358684068 -1.7758915240
C -1.8296740847 -1.1465175506 -0.5444949898
C -1.4268441422 -2.5633638543 -0.2308997241
C -1.6530482167 -0.1460785668 0.4082593299
C -2.0439493279 1.1630869836 0.1305154182
O -1.9305596711 2.1983235634 0.9990828078
H -1.2003629359 3.0217346306 2.7262971868
H -1.3943747537 1.2462085997 2.8296833104
H 0.7319447437 2.6783949005 1.0785449878
H 0.8086725407 0.1766513728 2.8640131241
H 3.2361043614 0.5779133975 2.5800608176
H 2.9429714280 1.6751678780 1.2168056076
H 3.7705997359 -0.5934415113 0.5290325792
H 2.3197285611 -1.3262388519 1.2303528490
H 2.3530861556 0.7287477483 -1.0429599793
H 0.8970543304 0.0117095184 -0.3400784601
H 1.6755236766 -2.2602423039 -1.0532884067
H 1.5119163852 -1.1570786327 -2.4287278953
H 3.1172183154 -1.5327563377 -1.7801746401
H -2.9237811573 2.4862797052 -1.3082386669
H -3.2380402552 0.7066147322 -3.0077839698
H -2.5339668703 -1.6104563479 -2.5222966220
H -1.0206620652 -3.0537871017 -1.1160991236
H -2.2983720908 -3.1326597682 0.1005183787
H -0.6800413687 -2.5911646697 0.5628480847
H -1.2196703602 -0.4019821983 1.3666930871
End
BondOrders
1 2 1.0
1 15 1.0
1 16 1.0
1 17 1.0
2 3 2.0
2 18 1.0
3 4 1.0
3 19 1.0
4 5 1.0
4 20 1.0
4 21 1.0
5 6 1.0
5 22 1.0
5 23 1.0
6 7 1.0
6 24 1.0
6 25 1.0
7 26 1.0
7 27 1.0
7 28 1.0
8 9 1.5
8 14 1.5
8 29 1.0
9 10 1.5
9 30 1.0
10 11 1.5
10 31 1.0
11 12 1.0
11 13 1.5
12 32 1.0
12 33 1.0
12 34 1.0
13 14 1.5
13 35 1.0
14 15 1.0
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
```